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Abstract 

A new technology was developed in this study which provides a successful numerical simulation of 
the whole process of flow transition in 3-D boundary layers, including linear growth, secondary insta- 
bility, breakdown, and transition at relatively low CPU cost. Most other spatial numerical simulations 
require high CPU cost and blow up at the stage of flow breakdown. A fourth-order finite difference 
scheme on stretched and staggered grids, a fully implicit time-marching technique, a semi-coarsening 
multigrid based on the so-called approximate line-box relaxation, and a buffer domain for the outflow 
boundary conditions were all used for high-order accuracy, good stability, and fast convergence. A new 
fine-coaise-fine grid mapping technique was developed to keep the code running after the laminar flow 
breaks down. The computational results are in good agreement with linear stability theory, secondary 
instability theory, and some experiments. The cost for a typical case with 162 x 34 x 34 grid is around 
2 CRAY-YMP CPU hours for 10 T-S periods. 


1 Introduction 

The transition process from laminar to turbulent flow in a wall-bounded shear flow is still a challenging 
and unsolved problem. Natural transition is a multi-stage process (Narasimha, 1990) involving 2-D linear 
evolution, 3-D secondary instability, breakdown, and transition (Figure 1). 

The linear stability equation was established by Orr (1907a, b) and Sommerfeld (1908), was solved by 
Tollmien (1931) and Schlichting (1932), and was experimentally confirmed by Schubauer and Skramstad 
(1948). The secondary instability was observed by Klebanoff, Tidstrom & Sargent (1962) for K-type 
fundamental transition, find was observed by Kachanov, Kozlov and Levchenko (1978) for subharmonic 
transition. The theoretical work was accomplished by Herbert (1983a & 1983b). There is really very 
little work, either theoretical or experimental, about the breakdown and transition zones which are the 
major parts of transition process. 


The numerical study of transition is still quite limited due to the lack of computational resources. 
First, most numerical studies are temporal (Orszag fc Kells, 1980; Wray & Hussaini, 1984; Kleiser & 
Laurien, 1985; Zang & Hussaini, 1986; Zang, Krist, Erlebader, & Hussaini, 1987). They can provide 
better resolution but lack physically realistic representation (Joslin et aL 1992). 

Second, although there have been some spatial studies (Fasel, 1976; Fasel & Bestek, 1986; Fasel & 
Konzelmann, 1990; Spalart, 1989; Danabasoglu, Biringen, & Streett, 1991), the spatial direct numerical 
simulation (DNS) is still in its early age (Kleiser fc Zang, 1991). Most of these can predict only the early 
stages of transition (pre-onset simulation) or fully developed turbulent flow without a transition process 
and require high CPU cost which is in the range of 100-1000 CRAY-YMP CPU hours. 


In contrast, the current study has two advantages. First, it was successful in spatial DNS for the 

whole process of transition including linear evolution, secondary instability, breakdown, and transition to 

turbulence. Second, the current approach is more efficient. The spatial DNS was carried out on a rather 

coarse grid (16 x 34 x 34 for each T-S wavelength) at an acceptable CPU cost which is in the range of 2 

- 10 CRAY-YMP hours. . . , 

boundary layer edge 



Figure 1. Idealized sketch of transition process on a flat plate. 
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2 Governing Equation In General Coordinates 


Let 


* — *(£»*?» C)» 
y = »(£,»?> C)> 
z = s(l,»bC)» 


the 3-D tim e-dependent incompressible Navier-Stokes equations can then be written as 

U = + v Zy + 

V = j(utj x + VI ly + wi?*), 

W=j(u£ e + vCy + wC z ), 


and 


du , r ,dUu dVu dWu 

lH + d£ + dr) + dC 

dv T (dU v , &V v . dWv 
dt dZ dr) d( 

dw , T/ dUw dVw dWw 

~di + J{ d£ + d 17 + dC 


) + ^ ^ 

) + (ty^ + “ ^ AlU 

) + (6^ + + &^) P " ^ Alty 

SU flF dW 
di + dr) + d( 


0, 

0, 

0 , 

0, 


( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 

( 6 ) 

(7) 


where it, v, w are velocity components, U, V , W are contravariant velocity components, P is pressure, Re is 
the Reynolds number based on the free stream velocity U^, the viscosity parameter v and some reference 
length, for example, $5 which is the displacement thickness of boundary layer at inflow, 


Re = 


U coSS 


and Ai is the physical Laplacian operator transferred to the computational (£, 77, C) space: 

i + yDiri + (Cl + + ©& + 2 (Z*v* + ZvVy + (,v*h 


d 3 


Ai = (£ + £y + + (vl + % + n z)g^ T IV. T v, T T -V*.*'/* T vy/y • sx-wy g^ 

Q 2 ft 3 

+ 2(£xCz + £yCy + + %{Vx(x + VyCy + *fcCx) + (£** + tyy + 

+ (l)xx + *Jyy + r) tI )— + (Cs* + Cyy + Cxx) • 


( 8 ) 


Here, we have 7 equations for 7 unknowns, u, v, u>, P, 17, V, and IF. 

The perturbation equations are obtained by decomposing the total flow into steady base flow and a 
perturbation. Using subscript 0 denote the base flow variables, and let 


v(x,y,z,t) *- vo (x,y,z)±tT(a?,y,z,i), 

V{x,y,z,t) <- V 0 (x,y,z) + V(x,y,z,t ), 

P(z,y,z,t) <- Po{x,y,z) + P(x,y,z,t), (9) 
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where v = (u,v,w), V — (U, V, W), and noting the base flow itself also satisfies the Navier-Stokes 
equations, we obtain the governing system for the perturbations: 


du , T/ d[u{u + U 0 ) + KoU] , d[u(v + V 0 ) + uoV] , 

, T , d[v(u + u 0 ) + v 0 u) , d\v{y + y 0 ) + v 0 v ] , 
at +l + drj 


a[t>(w + w 0 ) + voir] 


0 


a 


sc ) + ^ y ^ + ,7y aT ? + Cv ac )p D - AlV 


1 

Re' 


= o, 


= 0, 


dw , T/ d[w(u + u 0 ) + w 0 u] , a[w(v + y 0 ) + w 0 v ] , 

dt + 1 ^ + dTl 

du dv dw 

dt + dTj + dC 


0, 

0, 


( 10 ) 

( 11 ) 

(12) 

(13) 


Combined with (l)-(3), this system also has 7 equations and 7 unknowns for the perturbations. 


We perform the solving process as follows: 

1. Perform the surface and grid generation processes to obtain the required Jacobian coefficients. 

2. Solve system (1) — (7) to obtain the base flow solution* For a flat plate, we use the Blasius s imilari ty 
solution for the base flow. 

3. Solve system (1) — (3) and (10) — (13) to obtain the perturbation solution based on the above base 
flow. 


3 Boundary Conditions 

Benney-Lin type disturbances are imposed at the inflow boundary in this study: 

euRe<rf{<t>uue~ ia * } 

e3d+Real{<t>u3d + e'W z - ut) } 

t2 d Real{<t> v2d e«- a ¥- ut) } 
e 3d+ Real{4> v3d+ e^~ a ^ z - ut) } 

^RealiKzd^^-^} 

(14) 

where u> is the real frequency of the disturbance, /3 is a real constant that represents the spanwise 
wavenumber, and a = clr + iaj is the streamwise complex wavenumber obtained from linear stability 
theory. <f> u , <f> v and <f> w are eigenfunctions for the Orr-Sommerfeld equation. 


u(0,y,z,t) = 
+ 
+ 

/ As . 
v( — 2~,y,z,t) = 

+ 

+ 

+ 
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A no slip boundary condition is applied at the solid wall. According to the linear stability theory, the 
disturbances vanish at infinity, so we obtain the boundary conditions at far field given by 

u(x,y -* oo ,z,t) = 0, 
v(x,y-> oo, z,t) = 0, 

w(x,y —* oo, z,t) = 0. (15) 

Also, no pressure condition is needed at the inflow, solid wall, or far field since a staggered grid is used. 

4 Outflow Boundary Treatment 

Outflow boundary conditions have been the focus of study for the spatial simulation of flow transition 
by many researchers. For simplicity, we only describe the idea for a 2-D flat plate. The technique for the 
3-D problem is the same. 
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Figure 2. Extended computational domain. 


Taking the advantage of the staggered grid, we can obtain a fairly effective approach. First, a buffer 
domain technique developed by Streett & Macaraeg (1989) is applied to our problem. Thus, a buffer 
domain is appended to the end of the original outflow boundary to smear all possible reflections from 
the buffered outflow boundary (see Figure 2). The problem is, in general, that the conventional buffer 
domain is too long (usually four to eight T-S wavelengths), which greatly increases computation cost. Our 
goal is to maintain the accuracy in the original computational domain, and to eliminate all the possible 
reflection waves in a very short buffer domain. To realize the above goal, the governing equations in the 
buffer domain should be parabolicized to allow only strictly outgoing waves. Thus, a first buffer function 
6(f) is introduced here and applied to the streamwise viscous terms: 


d 2 U 

de 


HO 


d 2 u 

3f 2 ’ 


d 2 v 

d? 


HO 


d 2 v 

3f 2 ‘ 


(16) 


6(f) is a monotonically decreasing function that changes from 1 to 0 so that the upstream effects of the 
streamwise viscous terms will gradually disappear in the buffer domain. The essential feature here is that 
all damping mode disturbances at the buffered outflow boundary become zero. To understand this, we 
rewrite the first equation of (16) as 


HO 


d 2 U 

3f 2 


d 2 U 


d[ 




(17) 


5 


Since 6(f) -+ 0 at the buffered outflow boundary, and accordingly -jL- -» oo, then, we can also consider 

that the buffered outflow boundary is compressed from f = oo by the function %/&(£). Now clearly, if the 
disturbances are stable (damping modes), they will vanish at f = oo. To treat unstable (growing) modes, 
we need a second buffer function 6jk(f) that reduces the Reynolds number in the buffer domain gradually 
to less than the critical (or subcritical) Reynolds number and makes all the perturbation modes become 
damping 


Re Re 


( 18 ) 


Thus, the new modified governing equations in the computational (f , tj) plane become: 


dU 1 d(UU_ + 2UpU) d goV + yVo + gV . 
dt + Vn d V Vn * 

b Rt „ d 3 U 1 d 2 ,U. , d .U„, dP 


dV 




d{U 0 V + UV 0 + UV) d(2V 0 V + VV ) 


biu,. d 3 v^id 2 v dV dP 

-Ri5 (fc + ^ + 


9*1 ‘ 

dU 


di) 

dV 


di + dr, 


= 0, 


= 0, 


= 0. 


( 19 ) 


( 20 ) 

( 21 ) 


The buffer functions are chosen as follows: 

^original < f < Uolal , 

1 0 £ ^ ^ ^original j 

C ( ^~ L r gin — + 1 LvriginaX <t< Uo tah / 99 x 

^buffer ( **) 

1 0 < £ ^ L original • 

It is clear that the first function decreases from one to zero very rapidly as one moves from the original 
outflow boundary to the buffered outflow boundary. The second function increasing from 1 to c + 1 is 
a quadratic function that is continuously differentiable at the original outflow boundary. Note that the 
total effect of these buffer functions is that, toward the buffered outflow boundary: 

• the mo mentum equations become increasingly convection dominated in the direction, while the 
equations generally become parabolic; and 

• the momentum equations become highly diffusion dominated in the tj— direction. 

This treatment makes the outgoing waves propagate forward without reflection in the ^-direction, and 
any os cill ation in the 77— direction will be effectively smeared in the buffer domain. 
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Figure 3. Buffered outflow boundary points. 


Finally, we need to specify the buffered outflow boundary conditions under these modified governing 
equations. The parabolic character of the above equations requires only two boundary conditions. As 
m pnti™^ before, we have the disturbances tending to zero at the buffered outflow boundary (which is 
actually located at £ = oo), so one condition is 

P = 0. (23) 

This is a very important condition since the elliptic character of pressure has not been modified in our 
new governing equation. Any improper condition for P may cause trouble. For the second condition, we 
use the traditional extrapolation method for V , i.e., 


d 2 V 

d? 


= o. 


(24) 


Though this condition may not be so accurate, accuracy of solutions in the buffer domain is not so 
important, and the main concern is that there must be no reflection wave traveling back to the original 
computational domain. 


Referring to Figure 3, condition (23) is imposed directly on P by defining 

Pc = 0. 


(25) 


This is an implicit condition. With it, the discrete continuity equation associated with Pc is then used 
to define U at the buffered boundary: 

U E = U C - Vn ~ — A£. (26) 

At] 

Condition (24) is imposed by determining V at ghost points just outside the boundary: 

Vs = 2 V c - V w . (27) 


Note that the above treatment is only suitable for the perturbation equations. 
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5 Simplification 

In this study, we still use rectangular but stretched grids obtained by a special but relatively simple 
mapping (Figure 4): 


* = £, 

y - y{v), 

z = C‘ (28) 


This yields 

J = Vyi 
£* = (* = !■> 

£y = Cz = *?x = *?z — Cx = Cy = 0. (29) 


For our numerical simulation, we choose the transformation function 


y(v) = 


ymax&V 

V max& 4“ Vmaxijlmax ~ *}) 


(30) 


where Vmax is the height of the computational domain in the physical coordinates y, w is the height 
of the computational domain in the computational coordinate 17 , and a is a constant which can be used 
to adjust the concentration of grid points. This yields an inverse map 


v(y) 


T 7maxy(° r 4~ j/max) 
ymax{<T + if) 


(31) 


We can then obtain 


V*i 


%v 


Wrr iaxytMi^(^ 4~ Umax ) 
\jtmax& 4" Vmax {Vmax *?)P 

Irimaxvjcr + y mog ) 

2 /max( & 4“ 2/)^ 


(32) 

(33) 



Eigure 4. ^-direction stretched grid. 

Under the above mapping, the governing equations can be simplified: 

du , d[u{u + U 0 ) + no u] a[tt(V + Vo) + t*o V] , d{uW + i*o W] , , dP 

^ ^ + dC 

dv . d[v{u + U 0 ) + v 0 ^] , d[v(V + Vo) + V 0 F] , d[vW + v 0 W) dP 
K+Vyi Fi + Trj + ^ )+r h — 


= °> 

(34) 


(35) 
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dw ,d[w(U + U 0 ) + wqU] 

: %( + 


g[u>(F + Vo) + U> 0 F] aK H-wo^l . , 9P 


dr) 


9C 


*) + 

«e 


3C 


— -Aitu = 
Re 1 

0, 

(36) 

r dW 

0, 

(37) 

i + ~dC ~ 

u = 

u 

> 

(38) 


7 h 

w = 

w 
—— > 

(39) 


Vy 

V = 


(40) 


where ti, v, w, P, U, V, W are all fluctuating parts of the corresponding variables, and no, t>o> w o> Uo, Vo, Wo represent 
the base flow (Blasius solution for the flat plate). The transferred Laplacian operator in the computational space 
is simplified as 


a J , a J a’ a 

Al “ d? + Vy dp + dC 3 +7hv dr) 


(41) 


6 Discretization 

We use a uniform staggered grid for our problem in the computational {i,i),Q space (Figure 5). Letting <f> 
denote a generic function, the second-order backward Euler difference in time direction can be written as 

H 3 ^+* - 4 ^ + ^-* , 42) 

dt ~ 2Af 

and the fourth-order central difference in space can be written as 


-pj + 2AQ + *Pj + AQ - - A£) + P£ - 2Afl 

12 A£ 

df ~P£ + 2AQ + ljPj + 2A$) - 30^(Q + 16fl£ - Aj) - pj - 2Aj) 
dp U) ~ 12A£ 2 ’ 

d<f>„ . 1 A , x -pi + 2 A i) + 27 Pi + A i) - 27^(0 + Pi - Ai) 
di^ + 2 A <) 24A£ • 


(43) 

(44) 

( 45 ) 
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Figure 5. Staggered grid structure in the computational (£,*?, C) space. 


In the computational (£, 77, C) space, the grids are uniform. Suppose it, v, w and U, V, W are defined in 
terms of a staggered grid in the computational space (see Figure 5). Here, the values of P are associated with its 
cell centers, u and V with centers of the cell surfaces parallel to the (17, C) plane, v and V with centers of the cell 
surfaces parallel to the ((,() plane, and tu and W with centers of the cell surfaces parallel to the (£,ij) plane. 


Second-order backward Euler differences are used in the time direction, and fourth-order central differences are 
used in space. We can write the discretized governing equations symbolically as follows (Figure 6): 


Aeeuee + Aeue + Awuw + Aww^ww + Abn^nn + Anus + 

Asv.s + Assess + Affv.ee + Apup + Ab^b + Abb^bb ~ 

Ac vc + Dww Pww + DwPw + DePe ~ DcPc — S u , 
Beb v ee + Be v e + Bw v w + Bwwvww + Bnn*>nn + Bnvn + 

Bsvs + Bss^ss + Bff^ff + B f vf + Bbvb + Bbb'obb — 

Bcvc + Bs$Pss + EsPs + EjfPtr ~ EcPc = 

Ceewee + Ceu>e + Cwv>w + Cwwvww + + C N w N + 

Csvjs + Css w ss + Cff w ff + Cpwp + Cbidb + Cbbwbb ~~ 

Ccwc + FbbPbb + FbPb + PpPp ~~ PcPc — S w , 
DU eeB ee + DUeUe + DU W Uw — DTJcUc + DVhhVnn + 
+DV n V n + DV S V S - DV C V C + DW FP W FF + 

DWpWp + DW b W b - DW C W C = S M . 


(46) 


(47) 


(48) 

(49) 


As an illustration of the notation we use, relevant symbols for the discrete ^-momentum equation are depicted 
in Figure 4. The coefficients and source term for the interior points of the discrete momentum equation (46) 
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associated with uc are given as follows: 


Abe 

Ab 

Aw 

Aww 

Ann 

An 
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Ass 
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Ap 

Ab 

Abb 

Ac 

De 

D w 

Dww 


ic 


-wk? + 3 k^‘ +2V °-- ) ’ 

nap -§?<»+«.>. 

nap +^ wr +^ w ). 

“ 12ReA(* ~ I^ (r ' WW + 

-afe + ^04- + ’W-5B5- 

4< *c _ + 2yc 

ZReArj 1 3 A^ ( n+ ^ + ZReAr,’ 

4ag i ^ 1 h c ty i y \ ^ yc 

ZRe&r? + 3 Ar} [ ' + ’ ZReArj' 

-OT-^( v " +y »“ )+ ® 
+ +w °"' 

J_ + _L ( J_ + ^l + J_ ) 

2At 2JZe ' A£ 3 A tj 2 AC 3 '’ 


24 A$’ 
P C = 


27 

24A£’ 


24 A£’ 

— 4«c + lie" 1 . / — u O.N.wln» + 8«0 nV* — Stios^i + UOSS^JI 

2 A i + VyC{ 12 Atj 

—UQFFWff -f ZuofW} - %UqbW\> + Uq BB^ bb\ 

+ 12AC ' 


(50) 


Here, superscripts n and tl — 1 axe used to indicate values at previous time steps, and superscript to + 1, which 
indicates the current time step, is dropped for convenience. Lower case subscripts denote the approximate values 
of the v and w at points where the associated values of u with capital subscript are located (Figure 6). Other 
symbols used in the above formulas are as follows: 

a = n y, 7 = rtyy. (51) 
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Pigme 6. Neighbor points for ^-momentum equation 
(U are at the same points as u and are not shown here). 

All function values that are required at other than the canonical locations are obtained by fourth-order inter- 
polations in the computational space. For example (see Figure 7), 

V e = (9(Vc + Vn + Vtrw + Vw) - (Vsww +V nnww +Vsb + Vnite))/32, ( 52 ) 
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Figure 7. Neighbor points fox fourth-order approximation for V e . 

The coefficients for the tj— and C~ momentum equations are defined in an analogous way, and the discrete 
continuity equation is developed simply by applying the fourth-order central differences to each term. 


7 Approximate Line-Box Relaxation (ALB) 

7.1 2-D uniform grids 

The basic approximate box relaxation (AB) approach is to relax by boxes instead of points. With a 2-D uniform 
grid (Figure 8) as an example, we first describe the basic idea behind AB. 
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The generic form of the equations associated with a box for the 2-D uniform grids can be written as 


A%ue + Awvw + A n un + A s us A c uc + ^ 

” Sue > 

(53) 

Pc ~ Pb 

AgUgg + A^-UC + AffVWE + AgUSE AffUg+ ^ 

= ^u a > 

(54) 

B%VE + ByyVW + BtfVN + B S Vs B C VC + ^ 

= ) 

(55) 

Be v ne + B w vxw + B n vnn 4* B s vc B c vn + ^ 

— 

(56) 

UE — Uc VN “ Vc 

- 0. 

(57) 

Ax Ay 


Here, the superscripts represent the point at which the discretization is centered. We proceed in the box-by-box 
process with a few global point Gauss-Seidel relaxation sweeps on the momentum equations, changing u and v 
and holding P fixed. This means that the four momentum equations (53)-(56) in the box phase are approximately 
satisfied. Now, procee din g by boxes in some order, we perform distributed relaxation of the form: 


vc 

4 “ 

VC -€l, 

UE 

<“ 

VE +€ 2 , 

VC 

♦- 

vc — 6\y 

VN 

— 

VN 

Pc 


Pc + AP, 


(58) 


where the corrections are chosen to satisfy the discrete continuity equation and four discrete momentum equations 
associated with the box. Note that the old values of u,v, and P approximately satisfy the associated momentum 
equations, so we obtain the following system for the corrections, ei, €3, 61 , £2, and A P: 


(Agej + Ag ei )- — = 0, 

A P 

(Ag€j + -A^r€l) — — 0, 

AP 

(Bg6 3 + BgM- — =0, 

(*£*, + 2#*)- ^ = °> 

*1 + €3 + £3 _ c 

— a r — I — — — 

Ax Ay 


( 59 ) 


where 


y _ ,U£ -U C , V N - VC\ 

* ni — l * t * /' 

x Ax Ay 
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(59a) and (59b) yield 


(59c) and (59d) yield 


and together we have 


Therefore, (59e) can be written as 


or 


The correction are thus given by 


fl -Ac — 

T 3 - a= 


_ a _ - Bjf 

T 3 - p - ’ Bg-Bf 


£1 _ _ -Bg + Ay 

6l~ y ~ AC + 61' Ax' 

A C + a 


+ « . 61 + ^ 


Ax 


Ay 


!L — s 


7(1+ . ( 1 + ?) 5 l _ _ 

At~ + Ay ~ Sn 


Si = 


» 


7(1+1) , (Hj) 
Ax ' Ay 


r _ 

* - 

Cl = • 7> 

a 


and 


AP = (Agd + A%e 2 ) • A*. 

To simplify this scheme, note that for the case At << 1 , we have 



Aq ~ A§, 

Bg~B%, 


A% » A% y 

» Aw, 


Be » Bjf } 

B% » B<-, 

SO 

a 

w 1, 


P 

w 1, 


7 

Ay 

Ax 

For most cases, we have 

A$ > A%, 

Ac > Aw, 


> B% % 

B% > B$ , 

so a } (3 and can be approximated in 

general by 

A %~1 


or ~ 

Ag 


P ~ 

Bg . 

«r~ ■ 


7 ~ 

Bg Ay 
A% ’ Ax' 


It is these approximations and update formulas in (62) that we use in AB. 


(60) 


(61) 


(62) 
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7.2 General coordinates 

Foi general coordinates, the discrete momentum equations can still be written in the same generic form as that 
for Cartesian coordinates, but the continuity equation is changed to 

£g zEzl + Ye — V c = o. ( 63 ) 

Ay 

The physical velocities u and v have the following relations with the contravanant velocities U and V. 

U = au + bv, 

V = cu + dv. ' 

This leads to the discrete continuity equation written as 

a_EUE + beVE ~ ac^c - bcvc ctfuu + djywjv — cc^c - dcvc _ n 

- — — T 7 

Ay 

where the superscript ~ represents a point that is not located at a canonical position and therefore requires 
interpolation. Ass umin g that 

(bsAvs — bcAvc) < (asAus — acAuc), ( 64 ) 

(c/fAitjf — ccAuc) < (djfAvif—dcAvc), ( 65 ) 

then the correspond correction equation can be approximated by 

age 2 + ace i , dx6 3 4- dc$i _ q 

aI 

Note that the defining relations for ej, , and 


and 

01 ~ (iy+«c)7 ■ + *>) ' 

A t " r Aij 

where aE,a c ,dn,dc correspond to the mapping coefficients between u,v and U, V. 

7.3 Approximate line-box relaxation (ALB) for 3-D problems 

AB usually works well for 2-D problems, but frequently fails to provide fast convergence for 3-D problems. The 
basic idea of ALB is to satisfy the continuity equation for all boxes lying on one line simultaneously. Figure 9 gives 
the distribution of corrections in the (£, y) plane for the ALB. This kind of relaxation is very useful when the grids 
are anisotropic. Assuming for simplicity that a = P = 1, then according to Figure 9, ALB solving the discrete 
system (46)— (49) can be described as follows: 

• Freezing P, U , V t W, v, and tu, perform line Ganss-Seidel relaxation on (46) over tbe entire computational 
domain to obtain a new t*. 

• Freezing P, 17, V \ W , ti, and u>, perform line Ganss-Seidel relaxation on (47) over the entire computational 
domain to obtain a new v. 

• Freezing P, U, V, W, u, and t>, perform line Gauss-Seidel relaxation on (48) over the entire computational 
domain to obtain a new w. 

• Use transformation (38)-(40) to obtain new U, V, W. 


Ay 

62 can be expressed in the same form as for Cartesian coordinates: 
€2 

t=* 

1- 
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• For all j = 2,3, • “,71; — 1 at once: change *° sa ^ ls fy tie 

associated continuity equations, then update P,;* so that tne new l/, V, W and P as well as the associated 
transferred ti, v, w satisfy the three momentum equations. 


u i- 1 e t - f 6 

1 

Pie* 

-► c 

+ ap 6 

) — 

+ ^5 — ^6 

► 

^i+l 6 k + f 6 

U i-± 5 k ^5 

Pm 

+ ap 5 

— 

+ *4“*5 

_^«+| ik + { 5 

U i-$ 4k-e<_ 

Pm 

+ AP 4 
o — 

■f — £4 

4 k + C < 

J 3 Jfc " f 3 

j 

Pi3k 

V { ' k 

+ AP S 

D — 

-}- £3 — £3 

3 fc + e 3 

Pi-§ n- f 2_ 

P<2 *< 

+ APj 

D — 

2 k + c 2 






Figure 9. Distribution of corrections in the (£,*?) plane. 


Since all of the u y v } w have been previously relaxed, and the If, V y W are updated, we assume that equations 
(46)-(48) hold exactly. Let e, 6 , cr and A P represent the corrections for U } V, W and P, respectively. Thus, for 
cube ijk (see Figure 9), the correction equations corresponding to (46)-(49) are 


( A !+| ik^i+i Jlk + A i- 1 jk^i-i ; Ai> j 

= o, 

(67) 

£{ to - «f+i) - Cl to-* - w - C f 

= 0 , 

(68) 

(C;K> +i + Clvj ^ - C' 4 a ^ 

(jwjj* ik + ik)*i + (^wgf+i + 

= 0, 

(69) 

+i>Ci k(«i - «#+o - k(*#-t - to 

= » 

(70) 

J = 2, 3, — 1, 




where the superscripts represent the point at which the discretization is centered. This system has 4(n; — 2) 
equations for 4 (n; — 2) variables. Unfortunately, coupling between the correction variables makes the problem 
somewhat complicated. To develop a simpler approximate system, define 
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Then, equation (70) can be written in terms of the unknowns 5j only: 



Then we obtain the tridiagonal system 


as 63 


62 


i 3 k 

C3 d3 63 


62 


Sfn i 3 k 

*, 


* 

= 

; 

Cny -2 2 &»j -2 


2 


Srn i k 

Cfij — 1 &nj — l . 


- *n,-l . 


. i k . 


(75) 


Thus, Sj, j = 2, 3, ■ • • , », - 1 can be determined very efficiently. The other velocity corrections are given by 


Cj = u xj Sj, 

<Tj = Wzjtj, 

3 = 2j 3, • • • , 7ij — 1. 


The U, V, and W are then updated on all cells in the *, k tj— line as follows: 


P is then updated via 


u i+l i* 

«- 

j* +C J» 


Ui-l i* 

«- 

J* _ C j» 


W H fc+i 

«- 

W ij jfc+i + <Tj, 




w ijk-i-<rj> 
j = 2, 3, • nj — 1, 

(76) 



V”_i t +6,-i - S jt 
3 = 3,4 1. 

(77) 


Pijk Pijk + APj, (78) 

3 — - 1 . 


8 Semi-Coarsening Multigrid 

For the large-scale algebraic system of the 3-D flows that must be solved at each time step, the usual relaxation 
methods by themselves are much too slow. To obtain optimal efficiency, we use a multigiid scheme based on AB 
and ALB described in the previous section. For simplicity of discussion, we consider only the two-grid case. 

We use a full approximation scheme (FAS) to accommodate nonlinearities. A two-level FAS algorithm for an 
equation of the form 


L h 4> h = f h 


(79) 
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may be described loosely as follows: 


i) relax on L h 4> h = f h , 

ii) solve L 2h <f> 2h = L 2h ll h <f> h + Il h {f h - L h <f > h ), 

iii) replace $ h + I% h (4? h - 

The notation we have introduced includes the difference operators L h and L 2h y the restriction operators I^ h (for 
the approximation) and I 2h (for the residual), and the interpolation operator Jjj*. 

A full-coarsening strategy is generally ineffective for problems that favor special coordinate directions (e.g., 
anisotropic problems). To overcome this limitation, we consider now a special combination of semi-coarsening and 
line-box relaxation. The basic idea is to use line-box relaxation in one direction (say the y-direction) and coarsening 
only in the other two directions (x- and z-directions). A two-level staggered grid projected in the (x, y)-plane and 
(x,z)-plane is given in Figure 10. 

The full weighting restriction is still used here for transferring the residual from fine to coarse grids. The stencils 
can be expressed as follows: 

n h (R«) : f i \ 1 1 , 

L 8 4 a J 

ih k {R*) : f i * 1 . 

L 4 4 J 

' I 1 

: | | 

. a 8 

% k (Rm) : T f ( 80 ) 

L 4 4 j 

These stencils can be explained geometrically as shown in Figures 11-12. 
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Figure 10. Two-level staggered grid structure for x — z direction semi-coarsening: 
(a) fine grid projection on (x,y) plane, (b) coarse grid projection on (x,y) plane, 
(c) fine and coarse grid projection on (s, z) plane. 



Figure 12. Full-weighting restriction for y-momentum and continuity equations. 

For the restriction of variables, bilinear interpolation is used. Its stencils are 

£*(«) : [.j], 

n k (v) : [ | | ] , 

: [\ 3 ] » 

%\ p ) : [ | j ] . ( 81 ) 

For semi-coarsening, the coarse to fine transfer operators are based on linear interpolation: 

/^(A-u) : [ i ] for Au i OI [ I I ] fol Au 2’ 


19 













9 Fine-Coarse-Fine Grid Mapping 

The spatial DNS usually meets difficulties after the flow goes into the breakdown stage when the shear layer 
is developed and the vortex breaks down to small scale vortices. The numerical simulation will thus have a huge 
energy burst, and the disturbance velocity will be amplified by tens or hundreds of times somewhere inside the 
flow field . The code then blows up. Apparently, it is not the physical case, but is largely caused by that the 
grid we used is not fine enough to resolve the small eddies which play the role to generate dissipations. To keep 
the numerical simulation going, we developed a fine-coarse-fine grid mapping technique. To explain this technique, 
let us see what happens for a 1-D problem. We do the fine to coarse grid restriction and the coarse to fine grid 
interpolation at each time step: 

*r — r2fr, f oM 
u c - h u f 

u f - hh u c 

Here, l£ h is a linear restriction and I% h is a linear interpolation. 

sin(27ix/2L) sin(2nx/4L) sin(27vx/L) 



K < 1 modes K = 1 modes 

Figure 16. Fine-coarse-fine mapping. 


Define L as a section which has five grid points on the fine grid, and (assume ux = 0) 


So \uf d \dz 


as the amplification factor. 


Assume we have different frequency modes, e.g., 

• fir 2 * x ^ i 1 1 1 

sm(K T x), K= l,2’4-8’"' 

(K > 1 is not visible on this grid), we can approximately get a by numerical integration fox different modes. 

[sin(ir£f)+|sin(i^L)]f 

2 sm(Kir) + sin(2ijTir) 

sin(^) +sin(ifT) +sin(|iirx) + i sin(2ffx)’ 
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when K = 1, ol — 0, when K — ► 0, a — ► 1. 



Figure 17. Amplification factor for different modes. 

Figure 17 clearly shows that this kind of grid mapping, fine to coarse restriction and coarse to fine interpolation, 
significantly damps the highest frequency ( K = 1), but has only a little effect on other modes ( K < 1), and has 
almost no effect on low frequency modes ( K — * 0). 

The computational resources are still quite limited for DNS even if we use today’s largest supercomputer. For 
certain grids, the highest frequency which can be well simulated is K = 1. This highest mode may generate higher 
frequencies which can not be simulated by current grids and may cause the computation to fail. This fine- coarse-fine 
grid mapping damps the K = 1 mode and protects other frequency modes. Of course, we do not want to eliminate 
the K = 1 mode, but to restrict its energy growth. The actual procedure of this technique is 

1 . u c = Il h uf, 

2. ti?™ = (l - 0)uf d + I31$ h u c . 

Here, we choose 


+ w 2 ) } 

which is proportional to the perturbation energy. Therefore, there is very little damping to K = 1 when the 
perturbation is very small. In this way, we successfully keep the code running to simulate the whole process of 
transition: linear evolution, secondary instability, breakdown, and transition. Note that the large eddies play a 
much more important role in flow transition than do small eddies which correspond to high frequency modes. 
We have to sacrifice these small eddies due to the lack of computer resources. But, the physics of transition and 
turbulence are still simulated quite well due to the accurate representation of lower frequency modes corresponding 
to large eddies. 

10 Computational Results 

10.1 Comparison with linear stability theory (LST) 

To verify the accuracy of our approach, we compare our results with the linear theory by assuming a parallel 
steady base flow and imposing a small disturbance at inflow. The base flow is now uo(x, y) — uo(xo> y), v 0 (x, y) = 0, 
where tio(xo,y) is the Blasius similarity solution at inflow. Since it is easy for finite difference schemes to simulate 
the damping (stable) modes, we just test the growing (unstable) modes. 

Let Rcq = 900, Ft = 86 (w = 0.0774), and /? = 0 (2-D case). The Orr-Sommerfeld solution provides an 
eigenvalue 

a = ocr + taj = 0.2229 — *0.00451, 

and the associated eigenfunctions and <£J which are dependent only on y (Figure 18). 
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Strearmvise velocity eigenfunction <f>* 



—QJS l : _ » 

o 5 10 15 20 

Wall— normal velocity eigenfunction <fT 



V 

Figure 18. Stream wise and wall-normal velocity eigenfunctions of 
an unstable mode (i?e$ = 900) in 2-D flat plate flow. 

In LST, the disturbances are assumed to be traveling waves. In the 2-D case, this yields 

u = €e~ ajX {<f>j t co3(aRX — - <f>j3in(ccRX — w^f)), 

* = ee“ a/Z (<f> v R cos(a R x - u> R t) - ffisin^RX - w*t)). (83) 

The inflow boundary velocities can then be obtained by setting z = Xq for u and x — zq — for v (xo is the 
x— coordinate at inflow). We assume that e is small. Figure 18 shows the complex eigenfunctions for both u and v 
components <f>J } <j>% and <f>}) in the physical (x, y)~ plane, which are normalized by setting maxty^} = 1. 

A moderate computational domain is selected. In this case, the plate length (measured from the inflow location 
x = Xq) is set to eleven Tollmien-Schlichting (T-S) wavelengths, while the buffer dom ain is an additional single 
wavelength, making the length of the total computational domain twelve T-S wavelengths. Here we choose = 
75 and use a total computational grid with 362 x 50 grid points. The grid is highly anisotropic near the solid wall 
(Ax > Ay). 

Since we use a fully-implicit scheme, the time step is restricted mainly by accuracy consideration. For our tests, 
we take 


At = 


2x 

320u 


0.2537. 


The convergence rate of the semi-coarsening multigrid method, which we generally found to be about 0.2 per 
V(2,2) cycle (relax twice on each grid level before descending and ascending), is much better than the performance 
of single-grid relaxation, as shown in Figure 19. 
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continuity equation 



multi grid 
relaxation 


work units 


0 200 -*00 600 
Figure 19. Convergence history at a fixed time for multigrid and single-grid relaxation. 

In order to compare the computational results with the linear theory solution, which is accurate for the parallel 
wall-bounded base flow with small disturbances, we first assume that the base flow is everywhere given by 


Uo{t,ri) = Uo{0,v), ^,»7) = 0, e = 10~ 4 , 

where, J7 o (0, v) is the Blasius similarity solution in the computational (£, rj) plane at £ = 0 and that the displacement 
of the boundary layer is a constant, < 5 * 0 - Then Re*(= Re o),a, and w do not change along the streamwise 
direction. The streamwise and wall-normal velocity components, tt and v, of the disturbance after 13 T-S periods 
(t=13T) are compared with the solutions obtained by the linear theory at a vertical position close to the solid wall 
(y* = 1.3137 for tt and y* = 1.2448 for r) in Figure 20. Excellent agreement in both amplitude and phase between 
the computational results of our fourth-order finite difference scheme and the solution obtained by LST is observed 
in the physical domain. 



Figure 20. Comparison of the numerical and theoretical velocity components near the solid wall (y* — 1.3137 for 
u and y* = 1.2448 for t>). Re* = 900, Fr = 86, parallel wall-bounded base flow assumption is used grids: 362 x 50 
(11 T-S wavelength physical domain + 1 wavelength buffer domain). 
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To compare our numerical results with those obtained by LST more precisely, we check the profiles of different 
Fourier modes. 

The profiles of the disturbance waves are obtained by using the Fourier transformation. Generally, we can 
expand any continuous function in the spectral space. Since 

ft+T 

fi(Jfe) = ii/ u {x,y,t)e ikut dt, (84) 

2x Ji 

where u(x,y,t) is the complex disturbance velocity and u(fc) are the Fourier coefficients corresponding to the 
frequency Jfcu>, we can then obtain 


|«(*)| = 2y/a(k) 2 + b{ky , k = 0, 1, 2, ■ • ■ 

(85) 

ft+T 

(86) 

a(Jfc) = — / Real{u{x,y,t)}cos(kut)dt, 

2x Jf 

ft+T 

(87) 

b(k }= — / Real{u(x,y,t)}sin(ku>t)dt. 

J t 


Only the fundamental wave (Jfe = 1) is obtained in this case. Figure 21 depicts the streamwise and wall-normal 
disturbance velocity profiles at x* = 439.2 and 608.4 (x* = x/6' Xa ), and shows the excellent agreement with those 
obtained by linear stability theory. 



Figure 21. Comparison of the numerical and LST velocity profiles at z* = 439.2 and 608.4. 

To check the high order accuracy, two different grids (362 x 50 and 182 x 26) are tested. For this kind of 
unsteady problem, it is quite difficult to show the grid convergence. We define the relative L 2 error norm for both 
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u and v as 


iiz.iia = 

mu 


E«2 ’ 



where tt, v denote the numerical solution, and ti e , v e denote the LST solution. The results given in Table I show 
that the accuracy of our scheme under the above measure is about 0(h 3 * 5 ). 


grids 

IIZ.Ha 

IIZ.Ha 

182 x 26 

0.3315 

0.3243 

362 x 50 

0.0393 

0.0319 


Table I. Relative L 2 error norm for u and v after 
13 T-S periods (calculated in physical domain). 

10*2 Secondary instability and transition 

Now we return to the real world. The base flow is now not parallel to the solid wall, but the Blasius similarity 
solution. 

The computational domain is restricted to 

x € [zo,®o-M**o]> 

y € [0 yVmox], 

, x 3x 1 

z 6 l W\'W\ l 

where l x is the number of T-S wavelengths in the computational domain, and Ao is the T-S wavelength at inflow 
(the T-S wavelength A varies when the base flow is non-parallel). 

A Benney-Lin (1960) type disturbance is imposed on the inflow: 
tT(0,y, z,t) = 

where <j> 2 d{y) and <f>$d±{y) correspond respectively to 2-D and 3-D eignsolutions of the Orr-Sommerfeld equation 
and the superscript (Jfe) denotes different velocity components. Following is a typical case we chose: 


Re* 0 = 900, Fr = 86 (w = 0.0774), 

(3 — 0 . 1 , Umax — 50 , 

a 2d = 0.2229 - 0.00451i, 

azd = 0.2169 — 0.00419t, 

€2 d — 0.03, = 0.01, 

l x — 10 with 2 T-S wavelengths for buffer domain 

xq = 303.9, x en d = 593.6. 


The grid we used here is 162 x 34 x 34 (including eight wavelengths and a two wavelength buffer domain). The 
time step is set to ^ of tlie wave P eiio<i - 

It takes around 9 CRAY-YMP CPU hours for the code to run 30 T-S periods. Figures 22 and 23 depict the 
contours of relative helicity at different times which clearly show the process of K-type fundamental transition, 
A-wave formation, the peak and valley splitting and vortex breakdown. It is found that the breakdown begins 
at the second peak when the A-wave is intensified to certain degree and the shear flow is developed. The vortex 
breakdown further contaminates the flow field which leads to a transition process. The patterns of helicity at 
t =: 30T is very s imilar to those t = 7 T, which suggests that the whole process of transition has been built up after 
t = 7T. It turns out that less than 2 CRAY-YMP hours are needed to simulate the whole process of transition for 
a 162 x 34 x 34 grid and 7 T-S periods. Figures 24 and 25 give the contour plots of total vorticity magnitude on 
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the y = 0.0477 (x,z)-plane and spanwise vorticity on the z = 31.4 (x,y)-plane at different times which clearly 
show the process of vortex breakdown and formation of multiple spikes. The appearance of random moving small 
vortexes after breakdown provides a clue that the flow no longer keeps its laminar status. These contours also show 
a qualitative agreement with the laboratory experiment conducted by Saiic et al. (1984). 

The x-component of perturbation velocity, u, at different streamwise positions, x = 305,421, and 537, but 
with same y-coordinate (y = 0.3667) is shown in Figure 26. Although the disturbance imposed at inflow is a sine 
function as shown at x = 305, the perturbation at other points is largely amplified and distorted. The perturbation 
velocity no longer keeps its sine function shape, but starts oscillating very fast, showing that high frequency modes 
have been induced. 

We also averaged u and v on z = 31.4 (x,y)-plane at different streamwise positions, x = 305,363,421,479, 
and 537 after the transition process was built up. The time-averaged u and v are given in Figures 27, 28, 29, and 
30. Figures 28 and 29 depict the differences in S profiles between the Blasius similarity solution and computational 
results at x — 421 and 537, which qualitatively agree with the experimental results given by Suder, O’Brien, and 
Reshotko. Figures 28 and 29 also show that the tZ-profile in the transition zone is sharper than those of Blasius 
solution. The wall stress t — is then larger than that of laminar flows. 

Figure 30 gives the v-profile. The v is always positive in a laminar boundary layer. But, our computational 
results show that v varies from positive to negative and then becomes positive again. This is a typical sign that 
the flow is experiencing transition. Figure 31 shows that the spectrum of perturbation u becomes wider as the flow 
moves downstream. We also tried another case : 


Re* 0 = 732, 
u?2 d = 0.0909, 

u>$d = 0.04545, 

(3 0.2418, Umax — 50, 

a 2d = 0.2490- 0.00351i, 

a 3d = 0.1103- 0.00650f, 
e 2d = 0.015, e$ d ± = 0.005, 

l x = 10 with 2 T-S wavelengths for buffer domain 

x 0 = 248.2, x tnd = 437.5. 


The grid we used here is still 162 x 34 x 34 (including eight wavelengths and a two wavelength buffer domain). 
This was set to correspond to a subharmonic transition, but we still got a k-type fundamental transition when we 
chose e 2d = 0.03 and €$d± = 0-01* However, when we changed to e 2d = 0.015 and €3 d ± = 0.005, a subharmonic 
transition was very clearly observed ( see figure 32). 

There is really a lack of reliable experimental data for transitioned flow, which can be used to judge the 
computational results. Also, we need finer grids to get better resolution for post-onset flow or turbulent flow. We 
have to sacrifice those small eddies now. However, the results apparently provide physically correct simulations for 
the whole process of flow transition. 
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Figure 22. Front view of the relative helicity obtained on a 162 x 34 x 34 gird at different times. 
Re* = 900, Fr = 86,/3 = 0.1, e 2 d = 0.03, = 0.01. Flow direction is from left to right. 
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Figure 24 • Contour plots of the total vorticity magnitude obtained on a 162 x 34 x 34 
gird at different times on the y jj = 0.1123 (x,z)— plane. Re* = 900, Ft — 8 6,/3 = 0.1, 
€ 2 d = 0.03, € 3 ^ = 0.01. Contour interval is 0.02, flow direction is from left to right. 
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Figure 25- Contours plots of spanwise vorticity on the Zq = 0 (i, y)-plan 
at different times. Re* = 900, Ft = 86, /9 = 0.1, e 2 d = 0.03, e 3< * = 0.01. 
Contour interval is 0.02, flow direction is from left to right. 
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Figure 26. Temporal distribution of the perturbation velocity u on y = 0.367, x = 304,421, and 537. 

= 900, FV = 86,/? = 0.1, = 0.03, €&* = 0.015. 
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Figure 32. Bird view of the relative helicity of subharmonic transition obtained on a 162 x 34 x 34 gird at 
different times. Re* = 732, (3 = 0.2418, = 0.015, = 0.005. Flow direction is from left to right. 
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11 Concluding Remarks 

• The folly implicit time-marching and fourth-order finite difference scheme on a stretched and staggered grid is 
accurate enough to simulate the pre-onset transitional flow with a relatively coarse grid. The computational 
results agree with linear stability theory, secondary instability theory and some experiments. 

• The sim ulation with relatively coarse grids still can provide qualitatively correct prediction for transitional 
flow. It shows that the large eddies play more important roles in the process of flow transition. 

• The spatial DNS with relatively coarse grids meets trouble at the flow breakdown stage since the grid is 
not fine enough to resolve the small eddies which play the role to generate dissipations. We then need to 
introduce some numerical dissipation. The new fine-coarse- fine grid mapping technique can keep the DNS 
code running to simulate the whole process of transition, including the linear growth, secondary instability, 
breakdown, and transition. 

• To get more accurate DNS simulations, we still need finer gnds, especially for post-onset flows. 
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